// STL
#include <cmath>                       // std::nexttoward()
//#include <numeric>                     // std::numeric_limits<>
#include <vector>                      // std::vector<>

// jScience
#include "jutility.hpp"                // sort_indices_desc()

// stats++
#include "statsxx/postprocess/ROC.hpp" // ROC_pt, calc_ROC_curve()


//
// DESC: Efficient algorithm to generate ROC points (R).
//
// INPUT:
//
//     L       : set of test examples
//     f(i)    : probabilistic classifier's estimate that example i is positive
//
// NOTE: On OUTPUT, R is ordered by descending f(i).
//
// NOTE: In order to assign scores to the ROC points, recall the meaning of score:
// NOTE: ... For a continuous random variable X and a threshold (score) T, TPR and FPR are defined as:
//
//     TPR(T) = int_T^inf dx f_1(x)
//     FPR(T) = int_T^inf dx f_0(x)
//
// NOTE: ... For the given data:
//
//     T at (1,1) --- it is just the minimum of f(i).
//     T at (0,0) --- it is an infinitesimal amount above the maximum of f(i).
//
// NOTE: This implementation (including notation) is that of Algorithm 1 in:
//
//     T. Fawcett "An introduction to ROC analysis" Pattern Recognition Letters 27, 861--874 (2006)
//
// NOTE: ROC curves are insensitive to class skew. There is therefore no need to equalize P & N (number of positive and negative instances).
//
inline std::vector<ROC_pt> calc_ROC_curve(
                                          const std::vector<int>    &L,
                                          const std::vector<double> &f
                                          )
{
    std::vector<ROC_pt> R;

    // SORT L AND f BY DECREASING f VALUE
    std::vector<int>    L_srtd;
    std::vector<double> f_srtd;

    for( auto idx : sort_indices_desc(f) )
    {
        L_srtd.push_back( L[idx] );
        f_srtd.push_back( f[idx] );
    }

    int FP = 0;
    int TP = 0;

    // COUNT THE NUMBER OF POSITIVE (P) AND NEGATIVE (N) VALUES IN L
    double P = static_cast<double>( std::count(L.begin(), L.end(), 1) );
    double N = static_cast<double>( L.size() - P );

    // NOTE: f_prev is used for the purposed discussed below.
    // NOTE: ... It is also assigned to the (0,0) point.
//    double f_prev = std::numeric_limits<double>::lowest(); // T. Fawcett method
    double f_prev = std::nexttoward(f_srtd[0], (f_srtd[0] + 1.));

    for( auto i = 0; i < L_srtd.size(); ++i )
    {
        // NOTE: The purpose of this check is discussed in T. Fawcett.
        // NOTE: ... We want to process all values of equal score, before adding a point to the ROC curve.
        if( f_srtd[i] != f_prev )
        {
            ROC_pt tmp;
            tmp.FPR   = static_cast<double>(FP)/N;
            tmp.TPR   = static_cast<double>(TP)/P;
            // NOTE: As discussed above, it is f_prev that is assigned as the score here.
            tmp.score = f_prev;

            R.push_back(tmp);

            f_prev = f_srtd[i];
        }

        if( L_srtd[i] == 1 )
        {
            ++TP;
        }
        else
        {
            ++FP;
        }
    }

    // NOTE: Because TP and FP go P and N at different times (obviously), it is possible (probably) at the end points to be (1,y) or (x,1) where x or y != 1.

    // PUSH BACK POSITION (1,1)
    ROC_pt tmp;
    tmp.FPR   = static_cast<double>(FP)/N;
    tmp.TPR   = static_cast<double>(TP)/P;
    tmp.score = f_srtd.back();

    R.push_back(tmp);

    return R;
}
